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Abstract 

The advantages of using higher order, interpolatory basis functions are examined in 
the analysis of transverse electric (TE) plane wave scattering by homogeneous, 
dielectric cylinders. A boundarv-element/finite-element (BEM/FEM) hybrid 
formulation is employed in which the interior dielectric region is modeled with the 
vector Helmholtz equation, and a radiation boundary condition is supplied by an 
Electric Field Integral Equation (EFIE). An efficient method of handling the 
singular self-term arising in the EFIE is presented. The iterative solution of the 
partially dense system of equations is obtained using the Quasi-Minimal Residual 
(QMR) algorithm with an Incomplete LU Threshold (ILUT) preconditioner. 
Numerical results are shown for the case of an incident wave impinging upon a 
square dielectric cylinder. The convergence of the solution is shown versus the 
number of unknowns as a function of the completeness order of the basis functions. 

Index Terms - BEM/FEM , Electromagnetic radiation and scattering, higher order, 
numerical analysis 


I. Introduction 


Two-dimensional, higher order basis functions have been shown to be advantageous both 
in FEM formulations and in integral equation formulations. Salazar-Palma, et. al., [1] 
examined higher order vector bases in various closed FEM problems. Hamilton, et.al., 

[2] examined higher order bases in TM scattering from a perfect electric conducting 
(PEC), circular cylinder. In each of these studies, the general observed trend was that 
increasing the completeness order of the basis functions improved the accuracy for a 
given number of unknowns, or reduced the number of unknowns required for a given 
accuracy. These benefits, however, are dependent upon the smoothness or regularity of 
the solution [1], Prior to the common use of edge-based basis functions, Gong and 
Glisson [3] used linear and quadratic nodal basis functions in a 2D hybrid FEM/Method 
of Moments formulation. Later, Peterson, et. al., [4] employed higher order (quadratic 
for the TM case and linear for the TE case) scalar and vector bases in 2D scattering from 
heterogeneous objects with mesh termination provided by an absorbing boundary 
condition. It was shown that the increase from piecewise constant bases to piecewise 
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linear bases significantly improved the accuracy of the normal component at media 
interfaces. 


The finite element method is often the preferred approach for modeling inhomogeneous 
or anisotropic regions. However, in unbounded regions, some method of terminating the 
finite-element mesh must be employed. Absorbing Boundary Conditions (ABC) and 
Perfectly Matched Layer (PML) absorbers have been used extensively to simulate an 
open boundary condition on the FEM mesh. While these methods preserve the sparsity 
of the system, they are not without drawbacks. In both, the mesh must extend away from 
the surface of the scattering object an acceptable distance, and some experimentation to 
determine this distance may be necessary. The use of the PML absorber results in poor 
conditioning of the resulting system [6], and material values for the component layers 
must be optimized over the band of interest. The ABC method involves higher order 
derivatives of the fields at the mesh termination. In addition, the ABC and PML are 
approximate boundary conditions. Since this work is intended for insertion into a 
production code (EIGER [7], [8]), we chose the BEM approach for mesh termination in 
order to minimize sources of error that might arise from a user’s incomplete knowledge 
of the characteristics or limitations of the termination model. 

In this paper, higher order basis functions are applied in the hybrid BEM/FEM 
formulation. The problem is formulated using the vector Helmholtz equation for the 
electric field in the homogenous, isotropic interior region. The fields on the exterior are 
represented in terms of equivalent electric and magnetic currents on the FEM region 
boundary, and a null field condition on the electric field just inside the boundary is used 
to obtain an Electric Field Integral Equation (EFIE) there [9]. It is shown that use of a 
generalized Gaussian quadrature is an accurate and efficient method for handling the 
singular self-term arising in the 2D EFIE. Numerical results are shown for a TE- 
polarized plane wave incident upon a square dielectric cylinder. Convergence of the 
equivalent electric and magnetic currents along the boundary is examined as a function of 
the completeness order of the basis functions. The Quasi-Minimal Residual (QMR) 
method with an Incomplete LU Threshold (ILUT) preconditioner is used to solve the 
resulting sparse system. The effects of the fill level and tolerance specified in the ILUT 
preconditioner are compared as the order of the basis function increases. 


II. Problem Formulation 

The problem addressed is TE plane wave scattering by a dielectric cylinder (Fig. 1). 

Each side of the cylinder is of length d = 0.58 A 0 , where A 0 is the wavelength in free 

space. The TE case was chosen since the operators involved have features similar to 
those encountered in the 3D case. The interior and exterior formulations are discussed in 
detail below, followed by a discussion on the enforcement of continuity of the tangential 
electric and magnetic fields at the interface. 
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Figure 1 . TE scattering from a square, dielectric cylinder. 


A. Interior Formulation 

The vector Helmholtz equation is used to model the electric field in the interior region, 

VxVxE - k 2 E=-jco/J J, (1) 

where it has been assumed that the material is homogeneous and isotropic throughout the 
domain. It is also assumed that there are no impressed current sources inside the 
dielectric region, so the right-hand side of (1) vanishes. Curl -conforming, interpolatory 
basis functions are used to represent the electric field in the interior region: 

E(p)«f;F v ^ ) (p), peS, (2) 

7=1 

where the superscript “p” denotes the completeness order of the basis, and N E denotes 
the total number of unknowns associated with the interior electric field. The 
completeness order is the highest order such that all independent polynomial terms of that 
order or less are represented. The so-called “Whitney forms”, or “edge bases”, are used 
as the 0 th order basis functions. Basis functions of completeness order p>0 are obtained 
by multiplying the 0 lh order edge basis function by a p ll '-order modified or shifted 
Silvester interpolatory polynomial [10]. 

Equation (2) is inserted in (1), and the weak form of (1) is then obtained by the Galerkin 
method, resulting in 
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( 3 ) 


7=1 V S 


£ v, J(vxii^ ) . (vxn< p) ))rfs' -k 2 ds' + 

5 

J(Vxn ( / ) )xft‘' ) .n^ = 0, i = l,...,N e , 


dS 


where Cl) p) is the i th testing function, and is the f h source basis function. The curl- 

curl term in (1) has been modified through the use of the divergence theorem, resulting in 
the first term of (3) and the introduction of a boundary integral term. This boundary 
integral term provides for coupling of the magnetic field to the exterior region through 
enforcement of the boundary condition on the tangential magnetic field. It can be 
thought of as providing the source for the FEM problem, as becomes more apparent when 
it is expressed in terms of the interior magnetic field at the boundary: 


n f r ' 

2^v\ {(Vxay^.fVxft^)^' -k 2 \£l ( p •£l\ p) dS' 

7=1 \S S 7 

= jcoju J(rixH) d£, i = l,...,N e 

dS 


( 4 ) 


Note that the basis function for the tangential magnetic field is identical to that used for 
the electric current in the exterior formulation. This is discussed in more detail below. 


B. Exterior Formulation 


The exterior region is modeled using the Electric Field Integral Equation (EFIE). Using 
the equivalence principle, electric and magnetic currents are placed around the boundary 
of the cylinder, and the total tangential electric field is set equal to zero just inside the 
boundary: 


M „ 

nx 
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jcoA +VO + — VxF +nxE,.„ f =0, ptaS 
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( 5 ) 


Here, E inc is the electric field associated with the incident plane wave in the absence of 

the scattering object, n is the normal to the contour, and the “ t ” implies that the 
bounding contour is approached from the inside. The jump discontinuity of the curl term 
caused by the magnetic current M has been removed and presented explicitly in the first 
term. Therefore, the electric field at the observation point p due to the curl term arises 
only from the magnetic current away from point p [1 1], Furthermore, all current, field, 
and potential dependence upon the position, p , on the cylinder has been suppressed for 
notational simplicity. Since the sources are radiating in a homogeneous region, the 
potentials can be expressed using the free-space Green’s function, G(p,p') , and (5) 
becomes 


* 


4 


( 6 ) 


M(p) 
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where dS 0 indicates that the integral does not contribute at p' = p , and the prime notation 
indicates source point coordinates. 

On the boundary segments, piecewise polynomial, interpolatory basis functions are used 
to represent the magnetic currents, which are z-directed in the TE formulation: 

M(p)*zf;F 7 .n^(p), peas, (7) 

7=1 


where N M is the total number of unknowns associated with the magnetic current, and the 
superscript “/?” is again the order of polynomial completeness. V. denotes the 
coefficients of the magnetic current, which constitute a subset of the coefficients, V j , of 

the electric field (2) as discussed below, n (/,) (p) is simply the product of the //’-order 
modified Sylvester polynomial that interpolates the / 1 ’ unknown and a pulse function that 
is unity when p is within the support of the segment, and zero otherwise. This choice of 
basis function allows for straightforward enforcement of tangential continuity across the 
boundary, dS , as discussed below. 

The electric currents are also represented by piecewise polynomial, interpolatory basis 
functions [10]: 

J(p)*X/ / .A< / ' ) (p), pedS (8) 

7=1 

with = C . Ay P) (p) , 

where C ; . is the unit vector along the / h segment, and A ( ; ' ;) (p) is the //’-order, 

interpolatory scalar function that interpolates the f h node in the set of Nj degrees of 
freedom associated with the electric currents. It should be noted that electric current 
bases for interpolation points that coincide with segment endpoints span the support of 
the contiguous segments in order to ensure continuity of current across the segments. 
Bases that interpolate interior points have only a single segment as their support. While 
this choice of basis function provides a differentiable function within the element and 
ensures current continuity between segments, it also implies that the lowest completeness 
order for the electric current is p = 1 . This is in contrast to the basis representations for 
the electric field and equivalent magnetic current, where the lowest order is p= 0. 
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As in the interior formulation, the Galerkin choice for testing functions is used to obtain 
the weak form of (6) by dot multiplying both sides by , yielding 


where an integration by parts and the divergence theorem are used to transfer derivatives 
from the scalar potential onto the testing function. The number of testing functions in (9) 
corresponds to the number of electric current unknowns, Nj . This is a result of the 
strong continuity enforcement in which the degrees of freedom associated with the 
magnetic current are shared with the electric field along the boundary of the interior FEM 

region. Thus, the set of unknowns corresponding to the magnetic field, | , is a subset 
of the total set of unknowns for the electric field, {V E } . 

C. Accurate Integration of Self-Term 

In order to realize the expected improvements in accuracy afforded by the higher order 
bases, it is necessary to accurately integrate all matrix entries; this requires that special 
consideration be given to the integration of singular self-terms. 

In two dimensions, the Green’s function in (9) is given by [1 1] 


where D =|p — p'j . As the source point, p' , approaches the observation point, p , the 
Hankel function may be approximated by [12]: 




( 9 ) 


{A ( n P >E inc (r)Jf, pe^, , 



( 10 ) 
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where In y is Euler’s constant (0.5772. . .). To accurately integrate the logarithmic 
singularity, the generalized Gaussian quadrature rule introduced by Ma et. al., [13] is 
used. This rule exactly integrates linear combinations of the functions 

1, ln(x),x,xln(x), x 1 , x 2 ln(x), ... x w_1 ln(x) (12) 

on the interval (0,1) with an A-point rule. To implement this quadrature rule, the source 
segment is simply divided into two subsegments on either side of the observation point 
and each is reparameterized with the parameter origin at the observation point where the 
Green’s function is singular. 
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Figure 2. Integration using generalized quadrature rule. 


D. Continuity Enforcement 

To uniquely determine the fields in the exterior and interior regions, continuity of the 
tangential fields at the boundary must be enforced in some sense. As a result of the 
interpolatory basis functions selected for the interior electric field and the exterior 
magnetic current, the nodes associated with each are collocated along the boundary, dS . 
Since it is the tangential component of the interior electric field that is interpolated at 
these boundary nodes, and the magnetic current is z -directed, assigning a single degree 
of freedom for each pair of collocated nodes for E and M enforces strong tangential 
continuity of the electric field across the boundary. Similarly, the same basis function 
and degree of freedom used to model the electric current on the exterior is also used to 
represent the tangential magnetic field in (4). 


E. Matrix Structure and Sparse Solver 

The linear system of equations consists of a total of Nj densely populated rows 
corresponding to the EFIE (9) and Ne sparsely populated rows corresponding to the 
Helmholtz equation (4). This system was solved by direct LU decomposition for the 
smaller problems considered. For larger problems, memory and speed considerations 
necessitate the use of an iterative method. Without any type of matrix normalization, the 
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hybrid system is very ill-conditioned. Therefore, column scaling was used prior to the 
iterative solution process [14]. The Quasi-Minimal Residual (QMR) iterative method 
was then employed with an Incomplete LU Threshold (ILUT) preconditioner [15]. This 
preconditioner approximates the complete LU factorization as 


LU & A = LU , (13) 

where A is the system matrix formed from (9) and (4), and L and U are obtained by 

applying a sparsity pattern to the standard LU factorization. For each row in L and U, 
the sparsity pattern may be affected by two parameters, a fill-in factor, filllimit, and a 
tolerance factor, TOL. The fill-in factor specifies the maximum number of entries, in 
addition to the number of non-zero row entries of the original matrix, to be allowed in the 
corresponding rows of L and U. Only the largest fill_limit entries in a row of L or U are 
retained, and these sets are denoted by S\ L) and S \ U) , respectively, for the i"' row. The 
final sparsity pattern is then determined dynamically for each row of L and U by 
disallowing any positions in S] L) and S] U) for which the corresponding value is less than 
a specified tolerance. This tolerance is computed as the product of the average 
magnitude of the row entries and the user specified tolerance factor, TOL. The tolerance 
criterion for the matrix element to be included in the sparsity pattern of L is summarized 
by 


k |>^x 
I A jyw 


I 


jeS) 


(A) 


(14) 


m= 1 


where £ (j is the (i,j) entry in the matrix L, TOL is the user-specified tolerance, and N ( ^ ] 
is the total number of non-zero elements in the lower part (i .e.,j < i ) of the z th row. A 
similar condition exists for entries of U. 


III. Numerical Results 


A. Eigenvalues of a Square Conducting Cylinder 


In order to verify the functionality of the finite-element formulation and code, the 
generalized eigenvalues, k in (1), were computed for the case of a square, cylindrical 
conducting cavity with an interior relative permittivity of e r = 1.0. For this case, the 
relative permittivity was set to unity. The theoretical convergence rate for the 
eigenvalues is [1] 


X, 



< Ch 2(p+l) , 


(15) 


where h is the extent of the element, A., is the numerically obtained estimate of the i th 

eigenvalue, is the i th exact eigenvalue,/? is the order of polynomial completeness, and 
C is a constant that depends on the element type, but is independent of h. Figure 3 shows 
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the convergence versus the number of segments for basis completeness orders ranging 
from one to five. The slopes for the theoretical rates of convergence are shown for 
comparison in the lower left comer where 5 = 2(p+l). The convergence rates compare 
favorably with the theoretical rates. The p =5 curve appears to converge faster than the 
theoretical rate; however, the two data points obtained for this case are probably 
insufficient to reveal the final convergence rate. 
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Figure 3. Convergence of lowest eigenvalue for a square cylinder with bases of 
completeness order p on triangular elements. 


B. Dielectric Cylinder 


In this section, results are presented for the case of TE plane wave scattering by a square, 
dielectric cylinder. The relative permittivity of the cylinder is e r =2.2, and the length of 
each side of the square is 0.58 . In Fig. 4, the equivalent electric current is plotted as a 

function of the normalized distance around the square (Fig. 1), where the normalized side 
length is unity. The basis order p is the completeness order for the electric current bases. 
In the case p= 5, the entire finite element region comprises just two triangles. It can be 
seen that the p = 5 case with just 85 unknowns compares very well with the p= 1 baseline 
case with 15,000 unknowns. 
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Figure 4. Equivalent electric current on a square, dielectric cylinder. 


Figure 5 is a similar plot of the equivalent magnetic current along the boundary. In this 
graph, p is the completeness order of the magnetic current bases. Note that both the p = 0 
and p = 4 solutions appear close to the baseline over much of the plot. However, the p = 0 
case has nodes nearer the vertices where the magnetic current is singular, and the error is 
significant at these points. 
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Figure 5. Equivalent magnetic current on a square, dielectric cylinder. 


Figure 6 shows the convergence of the equivalent electric current as a function of the 
completeness order of the electric current bases. The result obtained from a p= 1 solution 
with 15,000 unknowns was used as the baseline. For comparison, curves showing the 
rates of convergence. 



where N is the number of unknowns, are inset in the lower left comer of Fig. 6. These 
rates correspond to the convergence of a closed domain, FEM solution with different 
orders of polynomial approximation, a , when the local errors are of the same order of 
magnitude [1], Since the elements in this study were of uniform length, and the error will 
most likely be greater near the cylinder comers, this latter condition is certainly violated. 
Nonetheless, there is a dramatic improvement in convergence of the p = 2 case compared 
to /?= 1 . The p = 3 case shows marginal improvement over p= 2, and orders higher than p = 3 
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do not appear to offer convergence benefits for the electric current in this problem. It 
should be noted that, due to the edge singularity in the magnetic current, the higher order 
bases are probably disadvantageous near the comers of the cylinder. This edge 
singularity likely causes a diminishing return in accuracy with higher basis orders for this 
problem. 

Figure 7 shows a similar convergence plot for the equivalent magnetic current as a 
function of the completeness order of the magnetic current bases. As in Fig. 6, curves 
showing rates of convergence that correspond to the same closed domain problem 
described above are inset in the lower left comer for comparison. As was the case with 
the electric current, the error in the magnetic current is almost certainly not uniform 
across the elements, but it can be seen that the slopes for the cases p— 1,2, and 3 agree 
well with the inset curves prior to leveling off. Furthermore, Fig. 7 reveals that there is a 
significant advantage associated with the use of higher order bases up to p = 3 or 4, even 
though the convergence rates are not sustained as the nodes approach comers. 

Comparing Figs. 6 and 7, there is a more pronounced improvement with basis order for 
the magnetic current than for the electric current. This may be due, in part, to the fact 
that the increase in unknowns for the trials associated with lower basis orders results in 
more nodes near the comer where the magnetic current is singular. However, there is a 
noticeable decrease in the slope of all curves for which p > 0, and this also is due to the 
inability of the smoother basis functions to model the singularity near the comers. This is 
revealed in Fig. 8, which shows a distribution of the error for two trials with p= 2. As the 
number of unknowns increase from 408 to 1488, it can be seen that the convergence near 
the comers is very poor. 
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Figure 6. Convergence of equivalent electric current, J , as a function of basis order. 
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Figure 7. Convergence of equivalent magnetic current, M, as a function of basis order. 
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Figure 8. Distribution of error as a function of u for two trials with p= 2. 


C. Iterative Solver 

For the numerical convergence results discussed above, all of the problems except the 
baseline could be readily solved via direct Gaussian elimination. However, various trials 
were run to investigate the suitability of the Quasi-Minimal Residual (QMR) iterative 
solver with an Incomplete LU Threshold (ILUT) preconditioner for the hybrid system 
defined by (4) and (9). The trials consisted of different values for the tolerance and 
maximum fill-in parameters for the ILUT preconditioner. In all trials, a column-norm 
scaling was implemented prior to the iterative solution process. 

It was found that there existed narrow ranges of the tolerance and fill-in parameters that 
permitted stable operation of the iterative solver with the resulting preconditioner. Table 
1 below shows the two input parameters, tolerance, TOL, and fill-in factor, fill_limit, as 
well as the number of nonzero entries in the resulting L and U matrices for a p = 0 case 
with 2,880 unknowns. The first row corresponds to a complete fill-in, and, hence, is 
equivalent to a sparse direct solve method. The second to last column shows the 
cumulative number of nonzero entries in L and U as a percentage of the total number of 
non-zero entries required for a sparse direct solution {fill limit =100% and TOL= 0). In 
the final column, the convergence results are stated as the number of two-term QMR 
iterations required for an absolute residual of less than 10' 6 . It should be noted that cases 
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in which the preconditioner failed were marked by excessively large norms of the initial 
residual after the first approximate inverse. 

It can be seen from Table 1 that, for the case of p=0 with 2880 unknowns, the resulting 
size of the preconditioner is approximately 75% of the complete fill-in required for a 
sparse direct solution. Table 2 contains results for a p = 3 case with 2496 unknowns. As 
expected, the p = 3 matrix is not as sparse as the p = 0 matrix. In addition, the resulting size 
of the preconditioner, approximately 85% of the full fill-in from a direct solver, is larger 
than that achieved for the p - 0 case. In [5] it was found that quadratic and linear bases 
resulted in similar fill-in for a nearly equal number of unknowns. In contrast, a marked 
increase in fill-in for the p = 3 case can be seen in trial 1 of Table 2 compared to trial 1 in 
Table 1. It should be noted that, in the cases presented here, no re-ordering of the system 
matrix was attempted. 

Two additional trials were run for the case p = 0 with 15,000 unknowns. These results are 
shown in Table 3, where it can be seen that the resulting ILUT preconditioner is about 
70% of the size of the complete, sparse LU factorization. 


Table 1 . Trials with p=0, N = 2880 unknowns. 


Trial 

TOL 

fill limit 
(%) 

L 

U 

(L+U) 

(%) 

Iterations 

1 

0 

100. 

306,870 

299,877 

100.0 

1 

2 

0 

10. 

169,661 

299,877 

77.4 

10 

3 

0 

5. 

157,061 

299,877 

75.3 

15 

4 

0 

4 

153,945 

292,408 

73.6 

22 

5 

0 

3 

147,951 

246,512 

65.0 

failed 

6 

0.0025 

100. 

296,486 

298,339 

98.0 

5 

7 

0.01 

100. 

270,195 

293,339 

92.9 

9 

8 

0.05 

100. 

168,767 

269,262 

72.2 

69 

9 

0.075 

100. 

131,146 

253,211 

63.3 

failed 

10 

0.01 

10. 

154,354 

293,339 

73.8 

17 

11 

0.02 

10. 

135,629 

286,760 

69.6 

failed 
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Table 2. Trials with p=3, N = 2496 unknowns. 


Trial 

TOL 

filljimit 

(%) 

L 

U 

(L+U) 

(%) 

Iterations 

1 

0 

100. 

503,396 

495,692 

100.0 

1 

2 

0 

15. 

388,212 

495,692 

88.5 

12 

3 

0 

10. 

354,744 

495,692 

85.1 

16 

4 

0 

7.5 

286,316 

453,016 

74.0 

failed 

5 

0.0025 

100. 

412,582 

486,011 

89.9 

16 

6 

0.005 

100. 

368,128 

478,775 

84.8 

37 

7 

0.0075 

100. 

309,825 

467,131 

77.8 

failed 

8 

0.001 

10. 

335,273 

491,409 

82.7 

31 

9 

0.05 

10. 

148,974 

435,175 

58.5 

failed 


Table 3. Trials with p=0, N = 15,000 unknowns. 


Trial 

TOL 

filljimit 

(%) 

L 

U 

(L+U) 

(%) 

Iterations 

1 

0 

100. 

3728130 

3689417 

100.0 

1 

2 

0.005 

10. 

1647754 

3551338 

70.1 

21 


III. Conclusions 


A FEM/BEM formulation using higher order, interpolatory vector basis functions was 
presented for application to TE scattering from dielectric cylinders. A similar 
formulation including results for scattering by 3D objects will be reported in the near 
future. Although the formulation was applied to a homogeneous cylinder, it can be 
readily adapted to inhomogeneous cylinders as well. A generalized Gaussian quadrature 
[13] was implemented and found to provide accurate and efficient integration of the 
singular self-terms. 
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In order to verify the FEM code, the convergence of the lowest eigenvalue of a square 
conducting cylinder was examined as a function of mesh density (//-refinements) and of 
basis completeness order (//-refinements). Numerical results were presented for the case 
of a TE plane wave impinging upon a dielectric square cylinder. The equivalent electric 
and magnetic boundary currents were plotted for the cases of low (p=0 for the magnetic 
currents) and high (p = 4 for the magnetic currents) basis completeness orders. Finally, h- 
p convergence was examined for the equivalent boundary currents with a 15,000 
unknown p = 0 (magnetic current) solution serving as a baseline. For both the electric and 
magnetic equivalent currents, the use of higher order bases improved the convergence 
significantly up to order p= 3. The convergence rates tended to subside, however, beyond 
a certain level of //-refinement. It is suspected that this is due to the edge singularity of 
the magnetic current. Near edges, achieving full advantage from the use of higher order 
basis functions will likely require equilibrated meshes [1] or the use of singular basis 
functions [16], 

An ILUT preconditioner, used in conjunction with QMR, was shown to provide solutions 
with very few iterations, but provided only modest reductions in the required memory 
when compared to a complete, sparse LU factorization. Adjusting the tolerance and fill- 
in parameters for the implemented ILUT preconditioner did not permit a reduction of fill- 
in at the expense of QMR iterations. Rather, there appeared to be a distinct threshold 
below which the preconditioner provided a very good approximation to the inverse of the 
original system, but above which the preconditioner became very poorly conditioned. 

The ILUTP method, in which the “P” stands for pivoting, is recommended in [14] as a 
remedy for unstable results with the ILUT preconditioner, but has not been applied in this 
study. 
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